#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <mpi.h>
#include <time.h>
#define	NDEBUG
#include <assert.h>

//#define __DEBUG
#define TEST_RANK			0
#define UNUSED_PROC			-1

#define RANDOM(min, max)	(rand() % ((max) - (min) + 1) + (min))

#define	INDEX_SRC_FNAME					1
#define	INDEX_DST_FNAME					2	

#define ROW_REQ_TO_E		0
#define ROW_REQ_FROM_W		1
#define ROW_REQ_TO_W		2
#define ROW_REQ_FROM_E		3
#define ROW_REQ_CNT			4

#define COL_REQ_TO_S		0
#define COL_REQ_FROM_N		1
#define COL_REQ_TO_N		2
#define COL_REQ_FROM_S		3
#define COL_REQ_CNT			4

#define VER_REQ_TO_S		0
#define VER_REQ_FROM_N		1
#define VER_REQ_TO_N		2
#define VER_REQ_FROM_S		3
#define VER_REQ_CNT			4

#define	VER_WN				0
#define VER_EN				1
#define VER_WS				2
#define VER_ES				3

#define NORTH				0
#define SOUTH				1
#define WEST				2
#define EAST				3

typedef struct {
	int coord_h; // the coordinates of the processor in the grid
	int coord_w;
	int dim_h;  // the dimension of the processor grid
	int dim_w;
	char *data; // the address of the grid data
	char *count; // the count of occupied neighbors
	int data_h; // the dimention of the data
	int data_w;
	int data_size;
	int grid_rank; 
	int rank[4]; // the rank of neighbors
	char *nbr_data[4];
	char *nbr_addr[4];
	char ver_to_n[2];
	char ver_to_s[2]; 
	char ver_from_n[2];
	char ver_from_s[2]; 
	char vertex[4];
	MPI_Request row_reqs[ROW_REQ_CNT];
	MPI_Request col_reqs[COL_REQ_CNT];
	MPI_Request ver_reqs[VER_REQ_CNT];
	MPI_Status row_status[ROW_REQ_CNT];
	MPI_Status col_status[COL_REQ_CNT]; 
	MPI_Status ver_status[VER_REQ_CNT];
	MPI_Comm comm;
	int tag;
} type_grid_info;

typedef struct {
	int h;
	int w;  
}  type_dim;

typedef struct {
	int h;
	int w;
	int iter;
} type_dim_iter;

int get_dim_size(int number, int *height, int *width)
{
	int i, diff, result;	
	int temp_diff;

	diff = number;
	for (i = 1; i <= number; i++) {
		if ((number % i) == 0) {
			result = number / i;
			temp_diff = (result >= i) ? (result - i) : (i - result);
			if (temp_diff < diff) {
				diff = temp_diff;		
				*height = (result >= i) ? i : result; 
				*width = (result >= i) ? result : i; 
			}	
		}	
	}
	return 0;
}

void print_grid_data(const char *data, int h, int w)
{
	int i, j;
	const char *p = data;

	for (i = 0; i < h; i++) {
		for (j = 0; j < w; j++) {
			printf("%d ", *(p++)); 
		}
		printf("\n");
	}
}

char *read_grid_data(const char *f_name, int *height, int *width, int *iteration, \
				type_dim *proc_dim,  type_dim_iter **data_size_vec)
{
	FILE *fp;
	int i, j, value, proc_size, proc_h, proc_w;
	int row, column, rank;
	char *data, *data_buf, *p;
	char **sub_addr;
	int *size_h, *size_w;
	type_dim_iter *size_vec;
	char *row_buf;
	int row_buf_size;
	int ave_h, rest_h, ave_w, rest_w;
	char *p_src;

    MPI_Comm_size( MPI_COMM_WORLD, &proc_size );
	//get grid dimention
	get_dim_size(proc_size, &proc_h, &proc_w);

	fp = fopen(f_name, "r");
	if (fp == NULL) {
		return NULL;
	}	
	// the first two numbers
	fscanf(fp, "%d", height);
	fscanf(fp, "%d", width);
	// the third number is the number of generations
	fscanf(fp, "%d", iteration);

	/* when problem size is less than processor size, then use less processors */
	if (proc_h > *height)
		proc_h = *height;
	if (proc_w > *width)
		proc_w = *width;

	proc_dim->h = proc_h;
	proc_dim->w = proc_w;
	proc_size = proc_h * proc_w;	
	*data_size_vec = (type_dim_iter *) malloc (sizeof(type_dim_iter) * proc_size);
	size_vec = *data_size_vec;

	// allocate memory after getting the height and width
	data = (char *)malloc(sizeof(char) * (*height) * (*width));	
	data_buf = (char *)malloc(sizeof(char) * (*height) * (*width));	
	sub_addr = (char **)malloc(sizeof(char *) * (proc_size));
	size_h = (int *)malloc(sizeof(int) * proc_h);
	size_w = (int *)malloc(sizeof(int) * proc_w);
	row_buf_size = sizeof(char) * (*width) * 3;
	row_buf = (char *)malloc(row_buf_size);

	// caculate the size for local matrix
	for (i = 0; i < proc_size; i++) {
		size_vec[i].iter = *iteration;

		row = i / proc_w;
		column = i % proc_w;
		size_vec[i].h = *height / proc_h;	
		size_vec[i].w = *width / proc_w;	
		if (row < *height % proc_h)
			size_vec[i].h += 1;	
		if (column < *width % proc_w)
			size_vec[i].w += 1;	

		size_h[row] = size_vec[i].h;
		size_w[column] = size_vec[i].w;
	}

#ifdef __DEBUG
	for (i = 0; i < proc_size; i++) {
		printf("rank %d: h %d, w %d\n", i, size_vec[i].h, size_vec[i].w);
	}
#endif

	// read grid data from file to data_buf
	p = data_buf;
	while (fgets(row_buf, row_buf_size, fp) != NULL) {
		i = 0;
		while (row_buf[i] != '\0') {
			if ((row_buf[i] == '0') || (row_buf[i] == '1')) {				
				*(p++) = row_buf[i] - '0';
			}
			i++;
		}
	}
	fclose(fp);

#ifdef __DEBUG
	printf("original data:\n");
	print_grid_data(data_buf, *height, *width);
#endif

	// copy sub matrix data of each processor from data_buf to data 
	ave_w = (*width) / proc_w;
	rest_w = (*width) % proc_w;
	ave_h = (*height) / proc_h;
	rest_h = (*height) % proc_h;

	p = data;
	for (i = 0; i < proc_size; i++) {
		row = i / proc_w;
		column = i % proc_w;
		
		// caculate the start address of each sub matrix in data_buf
		sub_addr[i] = data_buf + column * (ave_w);	
		if (column < rest_w) {
			sub_addr[i] += column;
		} else {
			sub_addr[i] += rest_w;
		}
		sub_addr[i] += row * (*width) * (ave_h);
		if (row < rest_h) {
			sub_addr[i] += row * (*width);
		} else {
			sub_addr[i] += (rest_h) * (*width); 
		}
		
		p_src = sub_addr[i];
		// copy each sub matrix from data_buf to data	
		for (j = 0; j < size_vec[i].h; j++) {
			memcpy(p, p_src, sizeof(char) * size_vec[i].w);
			p_src += *width;
			p += size_vec[i].w;
		}
	}

	free(data_buf);
	free(sub_addr);
	free(size_h);
	free(size_w);
	free(row_buf);
	return data;
}

// save grid data to file
void save_grid_data(const char *f_name, const char *data, int h, int w, \
			type_dim proc_dim, type_dim_iter *size_vec)
{
	FILE *fp;
	int i, j, row, column, proc_h, proc_w, proc_size;
	char *data_buf;
	const char *p;
	char **sub_addr;
	char *txt_buf, *p_txt;
	int txt_cnt;

	//get grid dimention
	proc_h = proc_dim.h; 
	proc_w = proc_dim.w; 
	proc_size = proc_h * proc_w;

	data_buf = (char *) malloc(sizeof(char) * h * w);
	sub_addr = (char **) malloc(sizeof(char *) * proc_size);
	txt_buf = (char *) malloc(sizeof(char) * h * (w * 2 + 1));

	// copy data from *data to *data_buf
	p = data;
	for (i = 0; i < proc_size; i++) {
		row = i / proc_w;
		column = i % proc_w;

		// caculate the start address of each sub matrix in data_buf
		sub_addr[i] = data_buf + column * (w / proc_w);	
		if (column < w % proc_w) {
			sub_addr[i] += column;
		} else {
			sub_addr[i] += (w) % proc_w;
		}
		sub_addr[i] += row * w * (h / proc_h);
		if (row < h % proc_h) {
			sub_addr[i] += row * w;
		} else {
			sub_addr[i] += (h % proc_h) * w; 
		}

		// copy each sub matrix from data_buf to data	
		for (j = 0; j < size_vec[i].h; j++) {
			memcpy(sub_addr[i] + j * w, p, sizeof(char) * size_vec[i].w);
			p += size_vec[i].w;
		}
	}
	

#ifdef __DEBUG
	printf("final results:\n");
	print_grid_data(data_buf, h, w);
#endif
	p = data_buf;
	p_txt = txt_buf;
	txt_cnt = 0;
	for (i = 0; i < h; i++) {
		for (j = 0; j < w; j++) {
			*(p_txt++) = *(p++) + '0';
			txt_cnt++;
			*(p_txt++) = ' ';
			txt_cnt++;
		}
		*(p_txt++) = '\n';
		txt_cnt++;
	}

	fp = fopen(f_name, "w");
	fwrite(txt_buf, sizeof(char), txt_cnt, fp);
	fclose(fp);

	free(data_buf);
	free(sub_addr);
	free(txt_buf);
}

// return the cell value for next generation
int get_next_val(int val_now, int occupied_cnt)
{
	int val_next;
	if(val_now == 1) { // the cell is occupied
		if ((occupied_cnt < 2) || (occupied_cnt > 3))
			val_next = 0;
		else
			val_next = 1;
							
	} else { // the cell is empty
		if (occupied_cnt == 3)
			val_next = 1;
		else
			val_next = 0;
	}
	return val_next;
}

// caculate the vertexes
int cacl_conway_vertex(char **count, char **data, int height, int width, \
									char **nbr_edge, char *vertex)
{ 
	int occup_cnt = 0;
	char *nbr_n, *nbr_s, *nbr_w, *nbr_e;

	nbr_n = nbr_edge[NORTH];
	nbr_s = nbr_edge[SOUTH];
	nbr_w = nbr_edge[WEST];
	nbr_e = nbr_edge[EAST];
	if ((height == 1) && (width == 1)) {
		occup_cnt = 0;
		occup_cnt += (nbr_n != NULL) ? nbr_n[0] : 0; 
		occup_cnt += (nbr_s != NULL) ? nbr_s[0] : 0; 
		occup_cnt += (nbr_w != NULL) ? nbr_w[0] : 0; 
		occup_cnt += (nbr_e != NULL) ? nbr_e[0] : 0; 
		occup_cnt += ((nbr_w != NULL) && (nbr_n != NULL)) ? vertex[VER_WN] : 0;
		occup_cnt += ((nbr_e != NULL) && (nbr_n != NULL)) ? vertex[VER_EN] : 0;
		occup_cnt += ((nbr_w != NULL) && (nbr_s != NULL)) ? vertex[VER_WS] : 0;
		occup_cnt += ((nbr_e != NULL) && (nbr_s != NULL)) ? vertex[VER_ES] : 0;
		count[0][0] = occup_cnt;

	} else if (height == 1) {
		/* cacl west one vertex */
		occup_cnt = data[0][1];
		occup_cnt += (nbr_w != NULL) ? nbr_w[0] : 0; 
		occup_cnt += (nbr_n != NULL) ? (nbr_n[0] + nbr_n[1]) : 0; 
		occup_cnt += (nbr_s != NULL) ? (nbr_s[0] + nbr_s[1]) : 0; 
		occup_cnt += ((nbr_w != NULL) && (nbr_n != NULL)) ? vertex[VER_WN] : 0;
		occup_cnt += ((nbr_w != NULL) && (nbr_s != NULL)) ? vertex[VER_WS] : 0;
		count[0][0] = occup_cnt;

		/* cacl east one vertex */
		occup_cnt = data[0][width - 2];
		occup_cnt += (nbr_e != NULL) ? nbr_e[0] : 0; 	
		occup_cnt += (nbr_n != NULL) ? (nbr_n[width - 1] + nbr_n[width - 2]) : 0; 
		occup_cnt += (nbr_s != NULL) ? (nbr_s[width - 1] + nbr_s[width - 2]) : 0; 
		occup_cnt += ((nbr_e != NULL) && (nbr_n != NULL)) ? vertex[VER_EN] : 0;
		occup_cnt += ((nbr_e != NULL) && (nbr_s != NULL)) ? vertex[VER_ES] : 0;
		count[0][width - 1] = occup_cnt;

	} else if (width == 1) {
		/* cacl north one vertex */
		occup_cnt = data[1][0];
		occup_cnt += (nbr_n != NULL) ? nbr_n[0] : 0; 
		occup_cnt += (nbr_w != NULL) ? (nbr_w[0] + nbr_w[1]) : 0; 
		occup_cnt += (nbr_e != NULL) ? (nbr_e[0] + nbr_e[1]) : 0; 
		occup_cnt += ((nbr_w != NULL) && (nbr_n != NULL)) ? vertex[VER_WN] : 0;
		occup_cnt += ((nbr_e != NULL) && (nbr_n != NULL)) ? vertex[VER_EN] : 0;
		count[0][0] = occup_cnt;

		/* cacl south one vertex */
		occup_cnt = data[height - 2][0];
		occup_cnt += (nbr_s != NULL) ? nbr_s[0] : 0; 
		occup_cnt += (nbr_w != NULL) ? (nbr_w[height - 1] + nbr_w[height - 2]) : 0; 
		occup_cnt += (nbr_e != NULL) ? (nbr_e[height - 1] + nbr_e[height - 2]) : 0; 
		occup_cnt += ((nbr_w != NULL) && (nbr_s != NULL)) ? vertex[VER_WS] : 0;
		occup_cnt += ((nbr_e != NULL) && (nbr_s != NULL)) ? vertex[VER_ES] : 0;
		count[height - 1][0] = occup_cnt;

	} else {
		// W N vertex
		occup_cnt = data[0][1] + data[1][0] + data[1][1];
		if (nbr_w != NULL)
			occup_cnt = occup_cnt + nbr_w[0] + nbr_w[1]; 
		if (nbr_n != NULL)
			occup_cnt = occup_cnt + nbr_n[0] + nbr_n[1]; 
		if ((nbr_w != NULL) && (nbr_n != NULL))
			occup_cnt += vertex[VER_WN];
		count[0][0] = occup_cnt;

		// E N vertex
		occup_cnt = data[0][width - 2] + data[1][width -1] + data[1][width - 2];
		if (nbr_n != NULL)
			occup_cnt = occup_cnt + nbr_n[width - 2] + nbr_n[width - 1]; 
		if (nbr_e != NULL)
			occup_cnt = occup_cnt + nbr_e[0] + nbr_e[1]; 
		if ((nbr_e != NULL) && (nbr_n != NULL))
			occup_cnt += vertex[VER_EN];
		count[0][width - 1] = occup_cnt;

		//  W S vertex
		occup_cnt = data[height - 2][0] + data[height -2][1] + data[height - 1][1];
		if (nbr_s != NULL)
			occup_cnt = occup_cnt + nbr_s[0] + nbr_s[1]; 
		if (nbr_w != NULL)
			occup_cnt = occup_cnt + nbr_w[height - 1] + nbr_w[height - 2]; 
		if ((nbr_w != NULL) && (nbr_s != NULL))
			occup_cnt += vertex[VER_WS];
		count[height - 1][0] = occup_cnt;

		// E S vertex
		occup_cnt = data[height -1][width - 2] + data[height - 2][width - 2] + \
								data[height - 2][width - 1];
		if (nbr_e != NULL)
			occup_cnt = occup_cnt + nbr_e[height - 2] + nbr_e[height - 1]; 
		if (nbr_s != NULL)
			occup_cnt = occup_cnt + nbr_s[width - 1] + nbr_s[width - 2]; 
		if ((nbr_e != NULL) && (nbr_s != NULL))
			occup_cnt += vertex[VER_ES];
		count[height - 1][width - 1] = occup_cnt;
	}

	return 0;
}

// caculate the edge data except the vertexes
// cnt is the width of the edge,
int cacl_conway_edge(int data_cnt, char *count, char *data, int data_displ, \
			char *nbr1, int nbr1_displ, char *nbr2, int nbr2_displ)
{
	int i, j;
	int occup_cnt, cnt1, cnt2, cnt;

	for (i = 1; i < data_cnt - 1; i++) {
		cnt = *(data + (i - 1) * data_displ) + *(data + (i + 1) * data_displ);
		cnt1 = 0;
		cnt2 = 0;
		for (j = i- 1; j <= i + 1; j++) {
			if (nbr1 != NULL)
				cnt1 += *(nbr1 + j * nbr1_displ);
			else
				cnt1 = 0;

			if (nbr2 != NULL)
				cnt2 += *(nbr2 + j * nbr2_displ);
			else
				cnt2 = 0;
		}
		occup_cnt = cnt + cnt1 + cnt2;
		*(count + i * data_displ) = occup_cnt; 
	} 
	
	return 0;
}

// update the grid data 
int update_grid_data(char *data, char *count, int len)
{
	int i, j;
	for (i = 0; i < len; i++) {
		data[i] = get_next_val(data[i], count[i]);
	}
	return 0;
}

// flag = 0, don't check the elements at the edge
// count the neighbors one time
int cacl_conway_cnt_once(char *count, char *data, int height, int width, int flag)
{
	int i, j, m, n;
	int occupied_cnt = 0;
	int n_occup_cnt, s_occup_cnt, w_occup_cnt, e_occup_cnt;
	char **ptr_data, **ptr_count;

	ptr_data = (char **)malloc(sizeof(char *) * height);
	ptr_count = (char **)malloc(sizeof(char *) * height);

	// initialize the pointer
	for (i = 0; i < height; i++) {
		ptr_data[i] = data + i * width;
		ptr_count[i] = count + i * width;
	}

	// caculate the grid data except the cells at the edge
	for (i = 1; i < height - 1; i++) {
		for (j = 1; j < width -1; j++) {

			//caculate the number of occupied elements in the neighbors	
			occupied_cnt = 0;
			occupied_cnt += ptr_data[i - 1][j - 1] + ptr_data[i - 1][j] + ptr_data[i - 1][j + 1];
			occupied_cnt += ptr_data[i][j - 1] + ptr_data[i][j + 1];
			occupied_cnt += ptr_data[i + 1][j - 1] + ptr_data[i + 1][j] + ptr_data[i + 1][j + 1];
			assert((occupied_cnt <= 8) && (occupied_cnt >= 0));

			// save the count	
			ptr_count[i][j] = occupied_cnt;
				
		}
	}

	if ( flag == 1) { // caculate the elements at edge
		// caculate north and south edge
		for (i = 1; i < width - 1; i++) {
			n_occup_cnt = 0;
			s_occup_cnt = 0;
			for (m = 0; m < 2; m++) {
				for (n = i- 1; n <= i + 1; n++) {
					n_occup_cnt += ptr_data[m][n];
					s_occup_cnt += ptr_data[height - 1 - m][n];
				}
			}
			n_occup_cnt -= ptr_data[0][i];
			s_occup_cnt -= ptr_data[height - 1][i];
			ptr_count[0][i] = n_occup_cnt;
			ptr_count[height - 1][i] = s_occup_cnt;
		} 

		// caculate west and east edge
		for (i = 1; i < height - 1; i++) {
			w_occup_cnt = 0;
			e_occup_cnt = 0;
			for (m = i - 1; m <= i + 1; m++) {
				for (n = 0; n < 2; n++) {
					w_occup_cnt += ptr_data[m][n];
					e_occup_cnt += ptr_data[m][width - 1 - n];
				}
			}
			w_occup_cnt -= ptr_data[i][0];
			e_occup_cnt -= ptr_data[i][width - 1];
			ptr_count[i][0] = w_occup_cnt;
			ptr_count[i][width - 1] = e_occup_cnt;
		} 

		// four vetexes
		occupied_cnt = ptr_data[0][1] + ptr_data[1][0] + ptr_data[1][1];
		ptr_count[0][0] = occupied_cnt;

		occupied_cnt = ptr_data[0][width - 2] + ptr_data[1][width - 1] + \
												ptr_data[1][width - 2];
		ptr_count[0][width - 1] = occupied_cnt;

		occupied_cnt = ptr_data[height - 2][0] + ptr_data[height - 1][1] + \
												ptr_data[height - 2][1];
		ptr_count[height - 1][0] = occupied_cnt;

		occupied_cnt = ptr_data[height - 1][width - 2] + \
						ptr_data[height - 2][width - 1] + \
						ptr_data[height - 2][width - 2];
		ptr_count[height - 1][width - 1] = occupied_cnt;

	} 


	free(ptr_data);
	free(ptr_count);
	return 0;
}

int set_persistent_comm(char **row_ptr, type_grid_info *grid_info, MPI_Datatype column_type)
{
	int i;
	int width, height, rank_n, rank_s, rank_w, rank_e;
	int flag;
	width = grid_info->data_w;
	height = grid_info->data_h;
	rank_n = (grid_info->rank)[NORTH];
	rank_s = (grid_info->rank)[SOUTH];
	rank_w = (grid_info->rank)[WEST];
	rank_e = (grid_info->rank)[EAST];
	

	// To East neighbor 
	MPI_Send_init(&row_ptr[0][width- 1], 1, column_type, rank_e, 
						grid_info->tag, grid_info->comm, &(grid_info->row_reqs[ROW_REQ_TO_E]));
	// From west neighbor
	MPI_Recv_init(grid_info->nbr_data[WEST], height, MPI_CHAR, rank_w,
						grid_info->tag, grid_info->comm, &(grid_info->row_reqs[ROW_REQ_FROM_W]));
	// TO West neighbor
	MPI_Send_init(&row_ptr[0][0], 1, column_type, rank_w, 
						grid_info->tag, grid_info->comm, &(grid_info->row_reqs[ROW_REQ_TO_W]));
	// From East neighbor
	MPI_Recv_init(grid_info->nbr_data[EAST], height, MPI_CHAR, rank_e,
						grid_info->tag, grid_info->comm, &(grid_info->row_reqs[ROW_REQ_FROM_E]));
	// TO North neighbor
	MPI_Send_init(&row_ptr[0][0], width, MPI_CHAR, rank_n, 
						grid_info->tag, grid_info->comm, &(grid_info->col_reqs[COL_REQ_TO_N]));
	// From south neighbor
	MPI_Recv_init(grid_info->nbr_data[SOUTH], width, MPI_CHAR, rank_s,
						grid_info->tag, grid_info->comm, &(grid_info->col_reqs[COL_REQ_FROM_S]));
	// TO South neighbor
	MPI_Send_init(&row_ptr[height - 1][0], width, MPI_CHAR, rank_s, 
						grid_info->tag, grid_info->comm, &(grid_info->col_reqs[COL_REQ_TO_S]));
	// From North neighbor
	MPI_Recv_init(grid_info->nbr_data[NORTH], width, MPI_CHAR, rank_n,
						grid_info->tag, grid_info->comm, &(grid_info->col_reqs[COL_REQ_FROM_N]));

	// send two vertexes to south neighbor
	MPI_Send_init(grid_info->ver_to_s, 2, MPI_CHAR, rank_s, 
						grid_info->tag, grid_info->comm, &(grid_info->ver_reqs[VER_REQ_TO_S]));
	// recv from north neighbor
	MPI_Recv_init(grid_info->ver_from_n, 2, MPI_CHAR, rank_n, 
						grid_info->tag, grid_info->comm, &(grid_info->ver_reqs[VER_REQ_FROM_N]));
	// send two vertexes to north neighbor
	MPI_Send_init(grid_info->ver_to_n, 2, MPI_CHAR, rank_n, 
						grid_info->tag, grid_info->comm, &(grid_info->ver_reqs[VER_REQ_TO_N]));
	// recv from south neighbor
	MPI_Recv_init(grid_info->ver_from_s, 2, MPI_CHAR, rank_s, 
						grid_info->tag, grid_info->comm, &(grid_info->ver_reqs[VER_REQ_FROM_S]));

	return 0;

}

int row_communication(char **row_ptr, type_grid_info *grid_info, MPI_Datatype column_type)
{
	MPI_Startall(ROW_REQ_CNT, grid_info->row_reqs);
	return 0;
}

int col_communication(char **row_ptr, type_grid_info *grid_info, MPI_Datatype column_type)
{

	MPI_Startall(COL_REQ_CNT, grid_info->col_reqs);

	return 0;
}

int ver_communication(char **row_ptr, type_grid_info *grid_info, MPI_Datatype column_type)
{
	// row communication first, then col communication	
	grid_info->ver_to_n[0] = grid_info->nbr_data[WEST][0];
	grid_info->ver_to_n[1] = grid_info->nbr_data[EAST][0];
	grid_info->ver_to_s[0] = grid_info->nbr_data[WEST][grid_info->data_h - 1];
	grid_info->ver_to_s[1] = grid_info->nbr_data[EAST][grid_info->data_h - 1];

	MPI_Startall(VER_REQ_CNT, grid_info->ver_reqs);

	return 0;

}

int parallel_conway(type_grid_info *grid_info, int iteration)
{
	int i, flag;
	char **row_data_ptr, **row_cnt_ptr;
	char *nbr2;
	int nbr2_displ;
	MPI_Datatype column_type; 

	row_data_ptr = (char **)malloc(sizeof(char *) * grid_info->data_h);
	row_cnt_ptr = (char **)malloc(sizeof(char *) * grid_info->data_h);

#ifdef __DEBUG
	if (grid_info->grid_rank == TEST_RANK) {
		printf("dim_h %d, dim_w %d, data_h %d, data_w %d\n", \
			grid_info->dim_h, grid_info->dim_w,grid_info->data_h,grid_info->data_w);
	}
#endif

	for (i = 0; i < grid_info->data_h; i++) {
		row_data_ptr[i] = grid_info->data + i * grid_info->data_w;
		row_cnt_ptr[i] = grid_info->count + i * grid_info->data_w;
	}

	// create new data type for west and east edge
	MPI_Type_vector(grid_info->data_h, 1, grid_info->data_w, MPI_CHAR, &column_type);
	MPI_Type_commit(&column_type);

	/* create persistent communication */
	set_persistent_comm(row_data_ptr, grid_info, column_type);

	for (i = 0; i < iteration; i++) {
		/* non-block communication between processors in a  column */
		col_communication(row_data_ptr, grid_info, column_type);
		/* non-block communication between processors in a row */
		row_communication(row_data_ptr, grid_info, column_type);
		
		/* caculate the grid data locally except the data on the edge */
		if ((grid_info->data_h >=3) && (grid_info->data_h >=3))
			cacl_conway_cnt_once(grid_info->count, grid_info->data, grid_info->data_h, grid_info->data_w, 0);

#ifdef __DEBUG 
	if (grid_info->grid_rank == TEST_RANK) { printf("\r\ncaculated middle data:\n");
			print_grid_data(grid_info->count, grid_info->data_h, grid_info->data_w);
		}
#endif

		/* Wait for row requests to complete */
		MPI_Waitall(ROW_REQ_CNT, grid_info->row_reqs, grid_info->row_status);

		/* caculate the west and east edge data */
		// caculate the west edge
		if (grid_info->data_w >= 2) {
			nbr2 = &row_data_ptr[0][1];
			nbr2_displ = grid_info->data_w;
		} else {
			nbr2 = grid_info->nbr_addr[EAST];
			nbr2_displ = 1;
		}
		if (grid_info->data_h >= 3) 
			cacl_conway_edge(grid_info->data_h, &row_cnt_ptr[0][0], &row_data_ptr[0][0], grid_info->data_w, \
						grid_info->nbr_addr[WEST], 1, nbr2, nbr2_displ);
#ifdef __DEBUG
		if (grid_info->grid_rank == TEST_RANK) {
			printf("\r\ncaculated west edge:\n");
			print_grid_data(grid_info->count, grid_info->data_h, grid_info->data_w);
		}
#endif

		// the east edge
		if (grid_info->data_w >= 2) {
			nbr2 = &row_data_ptr[0][grid_info->data_w -2];
			nbr2_displ = grid_info->data_w;
		} else {
			nbr2 = grid_info->nbr_addr[WEST];
			nbr2_displ = 1;
		}
		// only have one colum, don't caculate east edge
		if ((grid_info->data_h >= 3) && (grid_info->data_w >= 2))
			cacl_conway_edge(grid_info->data_h, &row_cnt_ptr[0][grid_info->data_w -1],
				&row_data_ptr[0][grid_info->data_w -1], grid_info->data_w, \
				grid_info->nbr_addr[EAST], 1, \
				nbr2, nbr2_displ);

#ifdef __DEBUG
		if (grid_info->grid_rank == TEST_RANK) {
			printf("\r\ncaculated east edge:\n");
			print_grid_data(grid_info->count, grid_info->data_h, grid_info->data_w);
		}
#endif

		/* non-block communication of vertexes between processors in a  column */
		ver_communication(row_data_ptr, grid_info, column_type);
		
		/* wait for column requests to complete */
		MPI_Waitall(COL_REQ_CNT, grid_info->col_reqs, grid_info->col_status);

		/* caculatethe north and south edge */
		// caculate the north edge
		if (grid_info->data_h >= 2) {
			nbr2 = row_data_ptr[1];
			nbr2_displ = 1;
		} else {
			nbr2 = grid_info->nbr_addr[SOUTH];
			nbr2_displ = 1;
		}
		if (grid_info->data_w >= 3) 
			cacl_conway_edge(grid_info->data_w, row_cnt_ptr[0], \
				row_data_ptr[0], 1, grid_info->nbr_addr[NORTH], 1, nbr2, nbr2_displ);

		// the south edge
		if (grid_info->data_h >= 2) {
			nbr2 = row_data_ptr[grid_info->data_h -2];
			nbr2_displ = 1;
		} else {
			nbr2 = grid_info->nbr_addr[NORTH];
			nbr2_displ = 1;
		}
		// only have one row, don't caculate south edge
		if ((grid_info->data_w >= 3) && (grid_info->data_h >= 2)) 
			cacl_conway_edge(grid_info->data_w, row_cnt_ptr[grid_info->data_h - 1], \
					row_data_ptr[grid_info->data_h - 1], 1, \
					grid_info->nbr_addr[SOUTH], 1, row_data_ptr[grid_info->data_h -2], 1);

#ifdef __DEBUG
		if (grid_info->grid_rank == TEST_RANK) {
			printf("\r\ncaculated north and south edge:\n");
			print_grid_data(grid_info->count, grid_info->data_h, grid_info->data_w);
		}
#endif
	
		/* wait for vertex requests to complete */
		MPI_Waitall(VER_REQ_CNT, grid_info->ver_reqs, grid_info->ver_status);
		grid_info->vertex[VER_WN] = grid_info->ver_from_n[0];
		grid_info->vertex[VER_EN] = grid_info->ver_from_n[1];
		grid_info->vertex[VER_WS] = grid_info->ver_from_s[0];
		grid_info->vertex[VER_ES] = grid_info->ver_from_s[1];

#ifdef __DEBUG
		if (grid_info->grid_rank == TEST_RANK) {
			printf("\r\nfour vertexes:\n");
			print_grid_data(grid_info->vertex, 1, 4);
		}
#endif
		/* caculate the four vertexes */
		cacl_conway_vertex(row_cnt_ptr, row_data_ptr, grid_info->data_h, grid_info->data_w, grid_info->nbr_addr, grid_info->vertex);

#ifdef __DEBUG
		if (grid_info->grid_rank == TEST_RANK) {
			printf("rank %d\n", grid_info->grid_rank);
			printf("north nbr:\n");
			print_grid_data(grid_info->nbr_data[NORTH], 1, grid_info->data_w);
			printf("south nbr:\n");
			print_grid_data(grid_info->nbr_data[SOUTH], 1, grid_info->data_w);
			printf("west nbr:\n");
			print_grid_data(grid_info->nbr_data[WEST], 1, grid_info->data_h);
			printf("east nbr:\n");
			print_grid_data(grid_info->nbr_data[EAST], 1, grid_info->data_h);

		}

		if (grid_info->grid_rank == TEST_RANK) {
			printf("\r\ncaculated vertexes:\n");
			print_grid_data(grid_info->count, grid_info->data_h, grid_info->data_w);
		}
#endif

		/* update grid data */
		update_grid_data(grid_info->data, grid_info->count, grid_info->data_size);
	}

	/* free the persistend requests */
	for (i = 0; i < ROW_REQ_CNT; i++) {
		 MPI_Request_free(&(grid_info->row_reqs[i]));
	}
	for (i = 0; i < COL_REQ_CNT; i++) {
		MPI_Request_free(&(grid_info->col_reqs[i]));
	}
	for (i = 0; i < VER_REQ_CNT; i++) {
		MPI_Request_free(&(grid_info->ver_reqs[i]));
	}

	MPI_Type_free(&column_type);
	free(row_data_ptr);
	free(row_cnt_ptr);
	return 0;
}

int set_mesh_topology(type_grid_info *grid_info, type_dim proc_dim, \
												int height, int width)
{
	//               not warp around,     reorder the ranks
	int dim_size[2], periods[2] = {1, 1}, reorder = 1;
	int orig_rank, grid_rank, orig_size;
	int dims[2], coord[2];
	int coord_h, coord_w, proc_h, proc_w;
	int dim_h, dim_w;

	MPI_Group orig_group, new_group;
	MPI_Comm new_comm;
	int *excl_rank, excl_size;
	int i, new_group_size;

    MPI_Comm_size(MPI_COMM_WORLD, &orig_size );
	MPI_Comm_rank(MPI_COMM_WORLD, &orig_rank);

	// get dimension of grid based on the processor size
	dim_size[0] = proc_dim.h;
	dim_size[1] = proc_dim.w;
	proc_h = proc_dim.h;
	proc_w = proc_dim.w;

	new_group_size = proc_h * proc_w;

	if (new_group_size < orig_size) {
		/* Extract the original group handle */
		MPI_Comm_group(MPI_COMM_WORLD, &orig_group);

		excl_size = orig_size - new_group_size;
		excl_rank = (int *) malloc(sizeof(int) * (excl_size));
		for (i = 0; i < excl_size; i++) 
			excl_rank[i] = i + new_group_size;
		MPI_Group_excl(orig_group, excl_size, excl_rank, &new_group); 
		MPI_Comm_create(MPI_COMM_WORLD, new_group, &new_comm);	

		if (new_comm == MPI_COMM_NULL) {
			free(excl_rank);
			return UNUSED_PROC;
		}
		MPI_Cart_create(new_comm, 2, dim_size, \
					periods, reorder, &(grid_info->comm));
		MPI_Group_free(&new_group);
		MPI_Comm_free(&new_comm);
		free(excl_rank);
	} else {
		MPI_Cart_create(MPI_COMM_WORLD, 2, dim_size, \
					periods, reorder, &(grid_info->comm));
	}
	
	MPI_Comm_rank(grid_info->comm, &grid_rank);
	MPI_Cart_get(grid_info->comm, 2, dims, periods, coord);
	
	coord_h = coord[0];
	coord_w = coord[1];
	dim_h = dims[0];
	dim_w = dims[1];
	grid_info->grid_rank = grid_rank;
	grid_info->coord_h = coord_h;
	grid_info->coord_w = coord_w;
	grid_info->dim_h = dim_h;
	grid_info->dim_w = dim_w;
#ifdef __DEBUG
	if (grid_info->grid_rank == TEST_RANK) {
		printf("new mesh topology:\n");
		printf("dim_h: %d, dim_w: %d\n", grid_info->dim_h, grid_info->dim_w);
	}
#endif
	grid_info->rank[NORTH] = ((coord_h + dim_h - 1) % dim_h) * dim_w + coord_w;
	grid_info->rank[SOUTH] = ((coord_h + 1) % dim_h) * dim_w + coord_w;
	grid_info->rank[WEST] = coord_h * dim_w + (coord_w + dim_w - 1) % dim_w;
	grid_info->rank[EAST] = coord_h * dim_w + (coord_w + 1) % dim_w;
	return 0;
}

int init_grid_info(type_grid_info *grid_info)
{
	grid_info->tag = 4242;
	grid_info->data = (char *)malloc(sizeof(char) * grid_info->data_size);
	grid_info->count = (char *)malloc(sizeof(char) * grid_info->data_size);

	grid_info->nbr_data[NORTH] = (char *)malloc(sizeof(char) * grid_info->data_w);
	grid_info->nbr_data[SOUTH] = (char *)malloc(sizeof(char) * grid_info->data_w);
	grid_info->nbr_data[WEST] = (char *)malloc(sizeof(char) * grid_info->data_h);
	grid_info->nbr_data[EAST] = (char *)malloc(sizeof(char) * grid_info->data_h);
	
	grid_info->nbr_addr[WEST] = (grid_info->coord_w == 0) ? \
							NULL : grid_info->nbr_data[WEST];
	grid_info->nbr_addr[EAST] = (grid_info->coord_w == grid_info->dim_w - 1) ? \
							NULL : grid_info->nbr_data[EAST];
	grid_info->nbr_addr[NORTH] = (grid_info->coord_h == 0) ? 
							NULL : grid_info->nbr_data[NORTH];
	grid_info->nbr_addr[SOUTH] = (grid_info->coord_h == grid_info->dim_h - 1) ? 
							NULL : grid_info->nbr_data[SOUTH];
	return 0;
}

int destruct_grid_info(type_grid_info *grid_info)
{
	free(grid_info->nbr_data[NORTH]);
	free(grid_info->nbr_data[SOUTH]);
	free(grid_info->nbr_data[WEST]);
	free(grid_info->nbr_data[EAST]);

	free(grid_info->data);
	free(grid_info->count);
	return 0;
}

int main(int argc, char *argv[])
{
	int i;	
	int err, orin_size, proc_size;
	double t_begin, t_end, t_par_begin, t_par_end;
	int w, h, iteration;	
	int dim_h, dim_w;
	char *grid_data, *grid_cnt;
	
	int orin_rank;
	int *send_cnt, *recv_cnt, *displs, recv_size;
	type_dim proc_dim;
	type_dim_iter *local_size_vec, local_size;
	type_grid_info grid_info;

    err = MPI_Init( &argc, &argv);
    err = MPI_Comm_rank( MPI_COMM_WORLD, &orin_rank );
    err = MPI_Comm_size( MPI_COMM_WORLD, &orin_size );

	/* rank 0 read data from specified file */
	if (orin_rank == 0) {
		if (argc != 3) { 
			// if you don't specify the name of the source and destination file
			printf("Usage: mpirun -np #proc <exe file> <data file> <result file>\n");
			return 1;
		} else { 
			t_begin = MPI_Wtime();
			printf("reading data from file...\n");
			/* read data from file */
			grid_data = read_grid_data(argv[INDEX_SRC_FNAME], &h, &w, &iteration, 
					&proc_dim, &local_size_vec);
			if (grid_data == NULL) {
				printf("Couldn't open %s, return\n", argv[INDEX_SRC_FNAME]);
				return 1;
			}
			proc_size = proc_dim.h * proc_dim.w;
			send_cnt = (int *)malloc(sizeof(int) * proc_size);
			recv_cnt = (int *)malloc(sizeof(int) * proc_size);
			displs = (int *)malloc(sizeof(int) * proc_size);
		}
	}	

	/* broadcast the processor dimention */
	MPI_Bcast(&proc_dim, sizeof(type_dim), MPI_CHAR, 0, MPI_COMM_WORLD);

	/* create a virtual grid topology */
	if (set_mesh_topology(&grid_info, proc_dim, h, w) == UNUSED_PROC) {
    	MPI_Finalize();
		return UNUSED_PROC;
	}

	/* scatter the grid data width and height to all processors */
	MPI_Scatter(local_size_vec, sizeof(type_dim_iter), MPI_CHAR, 
				&local_size, sizeof(type_dim_iter), MPI_CHAR, 0, grid_info.comm);	

	/* init grid_info data struct */
	grid_info.data_h = local_size.h; 
	grid_info.data_w = local_size.w; 
	grid_info.data_size = grid_info.data_h * grid_info.data_w;
	init_grid_info(&grid_info);

	/* Scatter grid data to each processor */
	if (orin_rank == 0) {
		send_cnt[0] = local_size_vec[0].h * local_size_vec[0].w;
		displs[0] = 0;
		for (i = 1; i < proc_size; i++) {
			send_cnt[i] = local_size_vec[i].h * local_size_vec[i].w;
			displs[i] = displs[i - 1] + send_cnt[i - 1];
		}
	}
	recv_size = grid_info.data_size;
	MPI_Scatterv(grid_data, send_cnt, displs, MPI_CHAR,
				 grid_info.data, recv_size, MPI_CHAR, 0, grid_info.comm);

#ifdef __DEBUG
	if (grid_info.grid_rank == TEST_RANK) {
		printf("scattered data of rank %d :\n", grid_info.grid_rank);	
		print_grid_data(grid_info.data, grid_info.data_h, grid_info.data_w);
	}
#endif
	
	/* parallel algorith */
	t_par_begin = MPI_Wtime();
	parallel_conway(&grid_info, local_size.iter);
	t_par_end = MPI_Wtime();

#ifdef __DEBUG
	if (grid_info.grid_rank == TEST_RANK) {
		printf("local results of rank %d :\n", grid_info.grid_rank);	
		print_grid_data(grid_info.data, grid_info.data_h, grid_info.data_w);
	}
#endif

	/* gather data from each processor */
	MPI_Gatherv(grid_info.data, grid_info.data_size, MPI_CHAR, 
			grid_data, send_cnt, displs, MPI_CHAR, 0, grid_info.comm);
#ifdef __DEBUG
	if (grid_info.grid_rank == TEST_RANK) {
		printf("final results sequential :\n");	
		print_grid_data(grid_data, h, w);
	}
#endif

	if (orin_rank == 0)	{
		/* save results to file */
		printf("saving results to file...\n");
		save_grid_data(argv[INDEX_DST_FNAME], grid_data, h, w, proc_dim, local_size_vec);		
		t_end = MPI_Wtime();

		// print the result
		printf("###Total run-time: %f s\n", (t_end - t_begin));
		printf("###Algorithm run-time: %f s\n", (t_par_end - t_par_begin));
		free(grid_data);
		free(send_cnt);
		free(recv_cnt);
		free(displs);
		free(local_size_vec);
	}
	/* destruct grid_info data struct */
	destruct_grid_info(&grid_info);

    err = MPI_Finalize();
    return err;
}
